function [equity_mat, debt_mat, lambda_mat, coupon_mat, optimal_firm, optimal_coupon, C_out] = equity_value_coupon(bParams,oParams,setts,upMat,f)
% Calculates the dynamic version of the equity value and also calculates rho

% Variables:
% bParams,oParams,setts		see readme
% upMat		update matrix equity value
% f			rate of creative destruction

% Set parameters
	pbar = oParams.pbar;
	lambda_max = setts.lambda_max;
    recovery_prec = setts.recovery_prec;

% No debt vs debt case
	if setts.nodebt==0
		C = setts.C;
	else
		C = 0;
	end
	Nc = length(C);

% Set initial values
	equity_mat = zeros(Nc,pbar+1);
	debt_mat = zeros(Nc,pbar+1);
	firm_mat = zeros(Nc,pbar+1);
	lambda_mat = zeros(Nc,pbar+1);
	coupon_mat = zeros(Nc,pbar+1);
	optimal_firm = zeros(1,pbar+1);
	optimal_debt = zeros(1,pbar+1);
	optimal_coupon = zeros(1,pbar+1);
    recovery_value = zeros(1,pbar+1);
    
% Itterate until recovery_value converges
    recovery_error = recovery_prec +1;
    while recovery_error>recovery_prec

% Save old recovery value
        recovery_value_old = recovery_value;
       
        % Reverse loop over coupon
        for i=Nc:-1:1


% Calculate new equity, debt, and firm value
            if i==Nc
                [equity_mat(i,:), debt_mat(i,:), firm_mat(i,:), lambda_mat(i,:), ~] = equity_value(bParams,oParams,setts,upMat,f,C(i),setts.equity_in,setts.debt_in,recovery_value);
            else
                [equity_mat(i,:), debt_mat(i,:), firm_mat(i,:), lambda_mat(i,:), ~] = equity_value(bParams,oParams,setts,upMat,f,C(i),equity_mat(i+1,:),debt_mat(i+1,:),recovery_value);
            end

% Save optimal coupon and firm value
            improvement = optimal_firm < firm_mat(i,:);
            optimal_coupon = optimal_coupon .* (1-improvement) + C(i) * improvement;
            optimal_debt = optimal_debt .* (1-improvement) + debt_mat(i,:) .* improvement;
            optimal_firm = max(optimal_firm,firm_mat(i,:));

% Output C
            C_out = C;
            
        end

% Recovery error
        recovery_value = optimal_firm;
        if setts.nodebt==1
            recovery_error = 0;
        else
            recovery_error = max(abs(recovery_value-recovery_value_old));
            disp(['Recovery value error is ',num2str(recovery_error,10)]);
        end                
	end

% Error optimal R&D lambda and C
    lambda_dum = lambda_mat(:,2:end,:);
    if setts.lambda_error==1 && max(lambda_dum(:)) == lambda_max
        disp('ERROR: Upper bound on lambda is reached');
    end

end
